Formation and evolution of density singularities in hydrodynamics of inelastic gases 
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We use ideal hydrodynamics to investigate clustering in a gas of inelastically colliding spheres. The 
hydrodynamic equations exhibit a finite-time density blowup, where the gas pressure remains finite. 
The density blowup signals formation of close-packed clusters. The blowup dynamics is universal 
and describable by exact analytic solutions continuable beyond the blowup time. These solutions 
show that dilute hydrodynamic equations yield a powerful effective description of a granular gas 
flow with close-packed clusters, described as finite-mass point-like singularities of the density. This 
description is similar in spirit to the description of shocks in ordinary ideal gas dynamics. 

PACS numbers: 45.70.Qj, 47.40.-x 



A central problem of non-equilibrium physics is to de- 
scribe clusters of matter developing from structureless 
initial states. Examples are numerous and include the 
large-scale structure of the Universe star formation 

, radiation-driven plasma condensations [|| , clustering 
of inertial particles in turbulent flows Q and clustering 
in inelastic gases 0, Q considered in this Rapid Commu- 
nication. 

Inelastic (or granular) g simple model of granu- 
lar flow Q, is a dilute assembly of hard spheres which 
lose energy at collisions. In the simplest version of the 
model the normal relative velocity of particles is reduced 
by a constant factor < r < 1 upon each collision. 
At collisions, internal degrees of freedom of the parti- 
cles absorb, apart from the energy, the entropy of the 
gas. As a result, inelastic gases exhibit a host of struc- 
ture forming instabilities, including the famous clustering 
instability of a freely cooling homogeneous inelastic gas 

1, IS 0, IE E3, EH \m EDI- This instability involves the 
development of macroscopic solenoidal flow (the shear 
mode) and potential flow (the clustering mode), the lat- 
ter causing the formation of clusters of particles where 
the particle density ultimately approaches the density of 
close packing of spheres The physical mechanism be- 
hind the initial density buildup is well understood [!, FLU ] . 
Collisional energy loss, enhanced by a local density ex- 
cess, can lead to a gas temperature decrease so strong 
that the local pressure falls, causing a further gas inflow 
and density growth. Pressure, therefore, plays a cen- 
tral role in the initial density build up, as corroborated 
by the linear theory of the clustering instability [!, @], 
where the density growth is accompanied by a decreasing 
in time pressure, having a local spatial minimum. How 
do non-linear effects change this linear scenario? Efrati 
et al. 12 1 and Meerson and Puglisi [3] addressed this 
question in the geometry of a long and narrow channel 
which we adopt in this work too. The channel geometry 
allows only one-dimensional (ID) macroscopic motions 
and therefore eliminates the shear mode. In the limit of 
strong instability (see below), Efrati et al. [T3| observed 
numerically a finite-time blow up of the gas density in 
ID hydrodynamic equations. Is it possible to establish 



the density blow up analytically, and does it develop for 
generic initial conditions? What is the role of pressure in 
the density blow up? Can hydrodynamics of dilute gas 
describe the state of the gas with already existing density 
spikes (interpreted here as close-packed clusters)? These 
questions are answered in this Rapid Communication. 

As clustering only intensifies with inelasticity, it suf- 
fices to consider the nearly elastic limit: 1 — r <C 1. This 
allows one (see Refs. [JEH) to employ hydrodynamics: 

d tP + d x (pv) = 0, p(d t v + vd x v) = -d x ( P T), (1) 
d t T + vd x T = -( 1 -l)Td x v-A P T 3 /\ (2) 

where p is the gas density, v is the velocity, T is the 
temperature (the particle mass is put to unity), p = pT 
is the pressure, and 7 is the adiabatic index (7 = 2 
and 5/3 in 2D and 3D, respectively). Equations (fT]) 
and ([2]) describe a ID dilute gas flow with bulk 
energy losses. These are accounted for by the term 
A/?T 3 / 2 coming from (1 — r 2 )T (proportional to the in- 
elastic energy loss per collision) times pT 1 / 2 (propor- 
tional to the collision rate). Kinetic theory yields A — 
2 7r (d-i)/2( 1 _ T ,2) .d-i /[dT(d/2)}, where d > 1 is the space 
dimension, and a is the particles diameter, see e.g. 0]. 
The omission of viscous and heat conduction terms in 
Eqs. ([T]) and ([2|) is justified in sufficiently long systems 
as, during its linear stage, the clustering instability pro- 
duces, at 1 — r <C 1, effective initial conditions with the 
Reynolds number Re ^> 1 

As in other ID compressible flows, it is convenient to 
go over from the Eulerian coordinate x to the Lagrangian 
mass coordinate m(x,t) = L p(x',t)dx'. Here m(x,t) 
measures the mass content between the origin [l9| and 
point x. This renders a simpler form of the equations: 



d t {l/p) = d m v, d t v = -d rn p, 
d t p = - 1P pd m v-Ap 3 / 2 p 1 / 2 . 



(3) 



Importantly, the cooling coefficient A can be elimi- 
nated from Eqs. ((3]) owing to the following scaling 
property: p(m,t) — p*(Am/A», At/A*), v(m,t) = 
D t (Am/A„ At/A*) and p(m,t) — p,(Am/A„ At/A*), 
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FIG. 1: The rescaled inverse density (a) and pressure (b) at 
t — 0.62 (circles), 0.77 (squares) and 0.99 (triangles). The 
solid lines show the analytical solution {5). (c): the inverse 
square root of the density at the singularity point m = 0. 



where it*, v* and solve Eqs. ([3]) for A* = 7v2. There- 
fore, in most of the following we put A = A* = jy2. 

Motivated by extensive numerical simulations (we used 
the classic artificial viscosity, staggered grid scheme of 
von Neumann and Richtmyer, see Ref. [18(1), we found 
a family of exact time-dependent analytical solutions of 
Eqs. ([3]), for which an initially smooth density profile 
blows up in a finite time t — t c . A distinct feature of 
these solutions is that the fluid elements have a constant 
in time acceleration, so v(m,t) depends on time linearly, 
while p(m, t) is time-independent. These solutions are: 



P*(m, t) 



p*(m, 0) 



[1 — ty/ Ap* (to, 0) cos m} 2 
(m,t) = — 2 



, p^ — 2Acosm, 



A cos m' . „ . . , , . 

- — -dm +2Atsmm, (4) 



/o V P*(m',0) 

where Pj (to, 0) is the initial density (an arbitrary func- 
tion) [l9(, and A > is an arbitrary constant. If p*(rn, 0) 
has a maximum at m = 0, then p*(m,t) blows up as 
(t c - t)' 2 at m = and t = t c = [Ap* (0, 0)}" 1 / 2 . As 
p must be non-negative, these solutions hold only at 
|m| < 7r/2, implying a finite mass of the gas m = tt. 
Depending on the behavior of p*(m,0) near |m| = 7r/2, 
this finite mass can be distributed over either an infinite, 
or a finite ir-interval. A simple example of the former, 

p*(m, t) — po costo(1 — tyj Apo cos m)~ 2 , p* = 2 A cos m, 
v*(m,t) = -2(A/p ) 1/2 m + 2 At sin m, (5) 

exhibits a density blow up at to = and t c — (Apo)~ 1//2 . 
We confirmed this solution numerically which indicates 
its stability with respect to small perturbations, see 

Fig. m 

Another exact solution (details will appear elsewhere 
jioj ]) shows that an emerging density blowup may coexist 
with an "ordinary" gas dynamic shock, and the Rankinc- 
Hugoniot conditions at the shock [l6| are obeyed. This 



solution involves, at t = 0, a zero-temperature gas 
with density p*(m, 0) = po exp (2^/j — 1 to) cos 3 m on 
the Lagrangian interval < to < 7r/2. The gas 
moves uniformly, v#(x, 0) = — i>o, and at t = hits a 
rigid wall located at x = 0. A reflected shock prop- 
agates into the gas (x > 0), leaving behind a flow of 
the form (J4| with p*(m,0) = p (-f + 1)[$ (to^cos to + 
$'(771)^/(7 - 1) costo] -2 and A = po«o(7 + l)/4- The 
function $(m) = J m exp [— m'ypy — l] dm'/ cos 2 to' de- 
termines the Lagrangian shock coordinate m s (t) implic- 
itly by t(m s ) = 2$(to s )/ [po^o(7 + !)]■ The collisional 
energy loss slows down the gas reflection from the wall, 
producing a density blowup at m = and t = t c = 

2V7^T/ [pov (7+l)j- 

These exact solutions (note that they are not self- 
similar) describe a density blowup ~ (t c — t)~ 2 at any 
A > 0. As p is finite (and non-zero), T vanishes at 
the singularity. The gas velocity is finite, while its x- 
derivative diverges as (t c — i) . In the Eulerian coor- 
dinate x the blowup of p occurs on a shrinking interval 
A(t) ~ (t c — t) 5 / 2 , while p* ~ x~ 4 / 5 and does not depend 
on time at A(i) -C \x\ -C 1. It is instructive to compare 
this new singularity with the well-known free-flow singu- 
larity [2]]], where the density blows up as (t c — t) -1 , the 
shrinking interval is A ~ (t c — t) 3 / 2 , and the power law 
tail of the density at singularity is ~ x~ 2 l z . For the free- 
flow singularity it is the Lagrangian velocity, rather than 
the Lagrangian acceleration, that is constant in time. 

Our numerical simulations with Eqs. |T]) and ([2|) 
showed that, for generic initial conditions, a finite-time 
density blow up always occurs [23|. Remarkably, the lo- 
cal flow structure near the singularity always coincides 
with that exhibited by the special solutions (j4|). One se- 
ries of simulations used the initial condition p(m, 0) = po, 
T(m, 0) = To and v(m, 0) = — vq tanh [to/ (pot)] at differ- 
ent values of two governing scaled parameters: the Mach 

1 /2 

number M = v /T ' and the cooling coefficient rescaled 
by PqI/M (and denoted by A with no ambiguity). The 
peak density p(0,t) exhibits a finite-time blow up, see 
Fig. [21 even at a relatively small A = 0.5, whereas p re- 
mains finite. Moreover, p _1 / 2 (0,i) goes to zero linearly 
as t — > t c . We observed these properties to hold at all A 
and M for which a high density is reached within a rea- 
sonable computation time (generally t c — > 00 as A — > 0). 
Figure shows that —9^ lnp at to = and t ~ t c vs. 
A is equal to A 2 /8: the prediction of Eq. (g]) with the 
restored A-dependence at 7 = 2. The same features 
(the density blowup and the universal flow structure near 
the singularity) are also observed when starting from a 
small-amplitude sinusoidal density or velocity perturba- 
tion around the homogeneous state. These results leave 
little doubt that the asymptotic behavior of the hydrody- 
namic fields near the singularity is described by Eqs. |4|. 

All the above implies that the non-linear evolution of 
the clustering instability, as described by the ideal hy- 
drodynamics of a dilute inelastic gas, brings about an 
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FIG. 2: The inverse density profiles at times 0, 2.06, 3.76, 
5.36, and 6.26 (a) and the evolution of the density (the solid 
line) and pressure (the dashed line) at the singularity point 
m = (b). The inset depicts p _1//2 (m = 0,t) close to t — t c . 
The parameters are y = 2, M = 1 and A = 0.5. 




FIG. 3: -d^ lnp(m, t) at m = and t ~ t c vs. A for M = 
(squares) and M = 10 (triangles). The solid line shows A 2 / 



infinite density. In sharp contrast with the linear stage 
of clustering instability [!, 0] , the gas pressure develops 
a local maximum at the point of the maximum density, 
and becomes time-independent in the Lagrangian frame. 
Obviously, to produce a density blowup, the flow must 
develop a region with d x v < 0. A negative shear at t = 
is necessary and sufficient for a blowup in the limit of 
A » 1, other parameters being fixed. This is because 
the pressure drops almost instantaneously there, result- 
ing in a free flow and a density singularity development 
in the regions where d x v < [rX[2l|. As the veloc- 
ity profile steepens, however, the compressional heating 
builds up and arrests the temperature drop, while the 
pressure again becomes important. The flow, having al- 
ready developed a high density peak, rearranges, and the 
blowup develops according to the new (non-free-flow) sce- 
nario. The rearrangement occurs within a time interval 
shrinking to zero as A — > oo. At the end of this inter- 
val the temperature dives to zero. So, as A — > oo, one 
has t c — > |min9 x v(;E,0)| , the blowup time of the free 
flow [2l[ . This explains the results of Ref . [l2j in a new 



light. Furthermore, if excluded volume effects are taken 
into account then, for sufficiently large A, the singularity 
remains of the free-flow type all the way until the close- 
packing density is reached [HI]. Still later dynamics is 



describable then by the Burgers equation [13, 23]. 

Now we return to our special solutions Q. Remark- 
ably, the dilute-gas equations (JTJ and © allow for an 
effective description of the states of inelastic gas with 
infinite density spikes. As differential equations of gas 
dynamics is indeterminate on a spike, one needs to re- 
turn to the integral form of the mass, momentum and 
energy balance. Let us continue the solution ((4]) beyond 
t = t c . We will describe this procedure for m > fl9j |. 
Assume that p*(m, t = 0) has a maximum at to = and 
is monotone decreasing at m > 0, vanishing at to = tt/2. 
Then at t > t c there is a point < m*(t) < Tt/2, implic- 
itly defined by 1 = ty/Ap^rnJ^), 0) cosm»(i), such that 

[m»(i),t], as prescribed by Eq. ([J}, is infinite. Now, at 
t < t c wc formally put m*{t) = so m*(t) is continuous 
at t = t c . Here is the solution in Eulerian coordinates, 
p(x,t), p(x,t) and v(x,t), which holds at any i>0: 

p(x,t) = 2m*(t)5(x) +p*(x,t), p(x,t) = p*(x,t), (6) 

where p*(x,t) and p*(x, t) are given implicitly at 
x > by the first line of Eqs. ([4]) and x(m,t) — 
Im (t) dm' / p*(m' , t). The gas velocity is 



Acosm' . . . . 

- — -am +2At [smm — smm,(f)] (7) 



m,(t)V p*(m',0) 



with the same x(m, t). We observe that, at t > t c , the 
solution includes a finite mass 2to*(£) concentrated at 
the origin. It is easy to see that Eqs. (JSj) and (JT]) solve 
Eqs. §3§ in the to— coordinate at to > m*(t) [but not at 
< to < m*(t)]. As a result, they solve Eqs. ((T|) and ^ 
in the x— coordinate at x > 0. It is left to verify that the 
integral form of Eqs. |T]) and © holds at x = 0. The 
mass flux pv, which is an odd function of x, has a dis- 
continuity at x = that produces the condition to* = 



lim 



p(x, i)v{x. t) 



lim 



m — >+m 



»(t) p{m,t)v(m,t) 



verifiable using L'Hopital's rule. The integral forms of 
the momentum and energy equations demand continuity 
of the momentum and energy fluxes at x = [for exam- 
ple, lim^o fl t p 2 T 3 / 2 dx = lim^oj^pp^^dx = by 
virtue of T(x = 0, t) = and finiteness of p(x = 0,i)]. 
It can be checked that Eqs. © and ([7]) obey these de- 
mands. As t — -> oo, all the gas mass accumulates at x = 0. 
The same continuation procedure, applied to our time- 
dependent solution with the shock, gives an example of a 
flow where a density spike and a shock are both present. 

According to molecular dynamic (MD) simulations 
[HI], the density buildup in the clusters is arrested only 
when the close-packing density is reached. An additional 
exact solution of Eqs. fl} and © establishes a simple 
relation between the density spikes and the close-packed 
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clusters. This solution involves a constant gas flux from 
infinity. In a moving frame the problem is equivalent to 
that of a piston entering, with speed vq, into a cold gas 
at rest. The latter problem was previously addressed by 



Goldshtein et al. |24( . Using an empiric equation of state 



where the pressure diverges at the close-packing density, 
they found that, at long times, the flow has three regions. 
First, there is a close-packed cluster, in contact with the 
piston; the cluster size grows linearly in time. It is fol- 
lowed by a region with non-trivial density, temperature 
and velocity profiles, separated from the gas at rest by 
a shock. The shock speed tends to vq in the dilute limit 
[24] ], Importantly, this flow is describable by the dilute 
equations JI]) and |2|). In the piston reference frame we 
obtain p = M(t)6(x)+p s t(x), v = v st {x) and T = T st (x), 
where p s t (x) , v st (x) and T st (x) form unique steady state 
solution of Eqs. {T]) and ([2]) obeying the boundary con- 
ditions p(x = +00, t) = po, v(x — +00, t) = —vq < and 
T(x = +00, t) — 0. In this description the close-packed 
cluster of Rcf. [1H becomes a single point x = which 
contributes the term M{t)8(x) to the gas density. The 
mass M(t) is coming from infinity, M = pqVo, and stored 
at x = 0. The shock is steady at x = L > 0, while at 
x > L the gas is at rest. In the region < x < L the 
mass and momentum fluxes are constant in space, so one 
can express v st (x) and T st (x) via p s t(x). The latter is 
given implicitly by 



x = - — < (7 + 3) arcsm v« , 

Apo I Vl - u 

where u = po/ Pst, and the A-dependence is restored. The 
piston-to-shock distance is 



L = 



1 

Apo 



1 



(7 + 3) arcsm ^ j-^ - ^8(7-!) 



This characteristic length scale behaves as A -1 , so it is 
very small at A > 1. The same feature holds for all 
solutions that we have presented here. 

In summary, ideal hydrodynamic equations for a freely 
cooling dilute inelastic gas describe the formation of 
density singularities in a finite time. What becomes 
of the singularities in microscopic theory? When the 
local gas density approaches the close-packing density, 
the singularities must get regularized and give way to 
(finite-density and finite-size) close-packed particle clus- 
ters [l3|, |T3| ■ One way to put the regularization into the 
continuum theory is to adopt an equation of state that 
diverges at the close-packing density [241] . An alterna- 
tive, suggested in this Rapid Communication, does not 
demand any regularization. It is similar in spirit to the 
ideal gas dynamic treatment of shock fronts (which in fact 
have a finite width) as discontinuities, that is point singu- 
larities of the first derivatives of the hydrodynamic fields 
[3]. As we have shown here, this alternative yields a 
powerful effective description in terms of smooth flow re- 
gions separated by point singularities: either "ordinary" 



shocks, or density spikes (close-packed clusters). It will 
be interesting to apply this description to driven flows of 
granular gases where theory can be compared with exper- 
iment. Finally, a direct comparison of our exact solutions 
with MD simulations is presently under way. 
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